function [worksheetName2Save,subfolder,helpDEBUG] = XXX_Function_ExecuteSaveMultiple_PETROGEN2020(infile,worksheetName2Save,subfolder)
tic
eval(infile)
%%% PROMPTED DECISIONS
%%%   Don't need to change anything here
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
WarnWave = [sin(1:.6:400), sin(1:.7:400), sin(1:.4:400)];
Audio = audioplayer(WarnWave, 22050);

play(Audio);
promptreponse = input(sprintf('\n Is ''%s'' the right folder to save the results? \n    Enter <strong>N</strong> to change, any key to continue: ',subfolder),'s');
if strcmp(promptreponse,'N')
    subfolder = input('Okay, what should it be?','s');
else
    disp('<strong>Okay! Using that folder</strong>')
    
end

disp('     ')

mkdir(subfolder);
addpath(subfolder)
warning('off','MATLAB:MKDIR:DirectoryExists')


play(Audio);
promptreponse = input(sprintf('\n Is <strong>''%s''</strong> the right .mat filename to save the results? \n    Enter <strong>N</strong> to change, any key to continue:',worksheetName2Save),'s');
if strcmp(promptreponse,'N')
    worksheetName2Save = input('Okay, what should it be?','s');
   
else
     disp('<strong>Okay! Using that folder</strong>')
end
disp('     ')
disp('<strong>Okay! Now running melting models...</strong>')





disp('      ')
disp('      ')



%% Don't worry about anything else below this line!

columnnum = 101;
addpath(genpath(pwd));
ColorsHERE


%%
% Import matrix with major element mantle bulk compositions
PetrogenMajorElement_Strings = {'SiO2','TiO2','Al2O3','Cr2O3','FeO','MgO','CaO','K2O','Na2O'};
DesiredDataName = 'USE';
[BulkCompositions, PetrogenBulkCompNames] = reorderElements_Rows(...
    'AAA_PETROGEN_SourcesAndDs', 'MantleBulkCompositions', DesiredDataName, PetrogenMajorElement_Strings,'FirstMajor','LastMajor');

Petrogen_BulkCompositionsNormalized = sumMgNumNormNoPT(BulkCompositions);


% Import matrix with trace element order and labels
[xlsNumbers, xlsText,xlsRAW] = xlsread('AAA_PETROGEN_SourcesAndDs', 'TraceElementOrder');
[ElementRow,ElementColumn]= find(strcmp(xlsRAW,'Elements')==1);
[LastElementRow,LastElementColumn]= find(strcmp(xlsRAW,'LastElement')==1);
PetrogenTraceElement_Strings = xlsRAW(ElementRow+1:LastElementRow-1, ElementColumn)';


% save indicies of important trace elements 
U_index = find(strcmp(PetrogenTraceElement_Strings,'U'));
Th_index = find(strcmp(PetrogenTraceElement_Strings,'Th'));
Ra_index = find(strcmp(PetrogenTraceElement_Strings,'Ra'));
H2O_index = find(strcmp(PetrogenTraceElement_Strings,'H2O'));

% Import trace element source compositions. Melt compositions can also be
% converted afterwards to have a difference source
[CoRow,CoColumn]= find(strcmp(xlsRAW,Petrogen_Co_Name)==1);
Petrogen_Co = xlsRAW(ElementRow+1:LastElementRow-1, CoColumn);
Petrogen_Co = cell2mat(Petrogen_Co);


% Import matrix with constant partition coefficients. Note we do change U
% and Th based on pressure later in the Petrogen code.
[DRow,DColumn]= find(strcmp(xlsRAW,'FirstD')==1);
[LastDRow,LastDColumn]= find(strcmp(xlsRAW,'LastD')==1);
Petrogen_Dmatrix = xlsRAW(ElementRow+1:LastElementRow-1, DColumn:LastDColumn);
Petrogen_Dmatrix = cell2mat(Petrogen_Dmatrix);

D_Th = eval(sprintf('Petrogen_Dmatrix(%d,:)',Th_index));
D_U = eval(sprintf('Petrogen_Dmatrix(%d,:)',U_index));
D_U./D_Th;





%%
runnumber = 1;
tic

MASTER_POOLS_INFO=[];
MASTER_POOLS_MAJORS =[];
MASTER_POOLS_TRACE =[];
MASTER_INSTANTS_MAJORS =[];
MASTER_INSTANTS_TRACE =[];
MASTER_INSTANTS_INFO =[];

MASTER_column_POOLS_TRACE =[];
MASTER_column_POOLS_MAJORS =[];
MASTER_column_POOLS_INFO =[];



MASTER_SM_MAJORS =[];
MASTER_SM_TRACE =[];
MASTER_SM_INFO =[];


%dont worry about major element residue for now:
%MASTER_Residue_MAJORS =[];
%MASTER_SolidResidue_MAJORS =[];
MASTER_Residue_TRACE =[];
MASTER_SolidResidue_TRACE =[];
MASTER_CPX =[];
MASTER_INSTANTS_Ds =[];
MASTER_POOLS_DISEQ=[];
MASTER_INSTANTS_DISEQ=[];
MASTER_SM_DISEQ = [];

RunConditions=[];

% set unused variables:
if strcmp(IsobaricOrPolybaric,'Isobaric')
    afewfunrnuns_U_vary=[NaN];
    TP_vary=[4000];
    afewfunrnuns_U_vary = [1];
    sU = {'NaN'};
    dFdP_Type = 'constant';
    GeodynamicModel='no';
    
    
    
    
else if strcmp(IsobaricOrPolybaric,'Polybaric')
        PIsobaric_vary=[NaN];
        sU = cellstr(strcat('U',num2str(afewfunrnuns_U_vary')))';
        sU = regexprep(sU,'U0.5','U05');
        sU = regexprep(sU,' ','');
    end
end


                        

%calculate the total number of model runs:
'Total Number of Runs'
totalNumberofRuns = ...
    size(dFdPwaterchange_vary,2)*...
    size(Porosity_vary,2)*...
    size(initialwatercontent_vary,2)*...
    size(TP_vary,2)*...
    size(afewfunrnuns_U_vary,2)*...
    size(BulkCompositions_vary,2)*...
    size(PIsobaric_vary,2)



for afewfunrnuns_dFdPwaterchange=1:size(dFdPwaterchange_vary,2)
    dFdPwaterchange = dFdPwaterchange_vary(afewfunrnuns_dFdPwaterchange);
    
    for afewfunrnuns_porosity = 1:size(Porosity_vary,2)
        porosity = Porosity_vary(afewfunrnuns_porosity);
        porosity_string = mat2str(porosity);
        porosity_string = regexprep(porosity_string,'\.*','');
        
        for afewfunrnuns_water = 1:size(initialwatercontent_vary,2)
            initialwatercontent=initialwatercontent_vary(afewfunrnuns_water);
            
            for afewfunrnuns_BulkCompositions = BulkCompositions_vary
                BulkCompositionHere = BulkCompositions(afewfunrnuns_BulkCompositions,:);
                BulkCompNameHere = PetrogenBulkCompNames(afewfunrnuns_BulkCompositions);
                
                
                % POLYBARIC LOOP
                for afewfunrnuns_temp = 1:size(TP_vary,2)
                    TP = TP_vary(afewfunrnuns_temp);
               
                    
                    
                    for afewfunrnuns_U = 1:size(afewfunrnuns_U_vary,2)
                        
                        
                        
                        if strcmp(GeodynamicModel,'yes')
                            clear Iv IT  Iv_2D Ivextend
                            infile = ['grid_',sU{afewfunrnuns_U},'_Tm',num2str(TP),'.mat'];
                            eval(['load ',infile])
                            % corrects the slope in the geodynamic model
                            
                            for kk = 1:size(IT,2)
                                IT_corrected(:,kk) =IT(:,kk)+ AdiabatSlope.*(Pb');
                            end
                            % extend geodyanmic model to whatever pressure you want
                            Pfin = 0.2;
                            Pinit = 30;
                            Pinc = 0.05;
                            Pb = Pinit:-Pinc:Pfin;
                            Pfin_extend = 30+Pinit;
                            Pb_extend = (Pfin_extend: -Pinc:(Pinit+Pinc));
                            Pbeek = [Pb_extend Pb];
                            Pb = Pbeek;
                            for kk = 1:size(IT_corrected,2)
                                Ttoextent(:,kk) =(max(IT_corrected(:,kk)).*ones(size(Pb_extend,1),1)+AdiabatSlope.*(Pb_extend-Pinit))';
                            end
                            Tb_extended_2D=[Ttoextent; IT_corrected];
                            
                            Tb_extended_onaxis = Tb_extended_2D(:,columnnum);
                            Tb = Tb_extended_onaxis;
                            
                            
                            Ivextend = repmat(Iv(1,:),size(Pb_extend,2),1);
                            Iv_2D = [Ivextend; Iv];
                            
                        end
                        
                        
                        if ismember(dFdP_Type, 'variable') == 1
                            VdFdP = ones(1,length(Pb))*0.002; %This will start off with a dFdP of 0.2%/kbar and increase to 2.5% at the surface
                            dFdP = VdFdP;
                        end
                        
                        
                        % ISOBARIC LOOP
                        for afewfunrnuns_Pisobaricchange = 1:size(PIsobaric_vary,2)
                            if strcmp(IsobaricOrPolybaric,'Isobaric')
                                
                                %set dFdP based on MaxF
                                Fcum = [0.01:.01:MaxF];
                                clear F Wnplus1S
                                Wo = 1;
                                Wnplus1S(1) =  Wo;
                                for uu = 1:size(Fcum,2)
                                    F(uu) = (Wo.*(Fcum(uu)-1)+Wnplus1S(uu))./Wnplus1S(uu);
                                    Wnplus1S(uu+1) = Wnplus1S(uu) - F(uu).*Wnplus1S(uu);
                                end
                                dFdP =F;
                                dFdP(isnan(dFdP)) =[];
                                
                                clear Pbadd
                                clear Pb
                                Pisobaricchange = PIsobaric_vary(afewfunrnuns_Pisobaricchange);
                                
                                Pb = Pisobaricchange*ones(size(dFdP));
                                Tb = TP*ones(size(dFdP));
                                
                                %This adds deeper pressures so that the mode is correct at the pressure
                                %of interest
                                Pbadd = 80:-0.1:Pisobaricchange;
                                Pbadd(end)=[];
                                dFdPadd = dFdP(1).*ones(1,size(Pbadd,2));
                                Tadd = zeros(1,size(Pbadd,2));
                                
                                Tb=[Tadd Tb];
                                Tb_extended_2D=repmat(Tb',1,201);
                                
                                dFdP=[dFdPadd dFdP];
                                Pb = [Pbadd Pb];
                                dummyPb = size(Pbadd,2);
                                
                                Ivextend=0.*Tb;
                                Iv_2D = 0.*Tb_extended_2D;
                                
                                dz=154.6073; %154.6073 m
                                dx=1; %1 m
                                x=-200000:2000:200000; %m
                                
                            end
                            
                            
                            
                            
                            runnumber
                            
                            
                            [XM_OUT,XOW_OUT,BAI3,BAI4,Ts,Tsdry,T,CumFrac,...
                                IncFrac,Fc,rXOW,MgRatio, NaKratio,test,initialBulkComp,...
                                CL_OUT,CS_OUT, CR_OUT,CL_BAI3,...
                                CS_minerals_OUT,Do_out,P_out,Dn_Test, ln, WnS, Ln, Qn, Lnr, Lnexp, Wnr,...
                                Gar2SpTransition_index, Sp2PlagTransition_index,Gar2SpTransition, Sp2PlagTransition,...
                                Lnexp_IT,DUTh_out, KDs,dFdPOUT,helpDEBUG] = ...
                                ...
                                ...
                                ...
                                PETROGEN2020(...
                                Tb,Pb,dFdP,...
                                BulkCompositionHere,...
                                porosity,...
                                dFdPwaterchange,initialwatercontent,...
                                Petrogen_Co,H2O_index, Petrogen_Dmatrix,U_index, Th_index,...
                                Tb_extended_2D,'initialTransition',forceSpPlag);
                            
                            
                            
                            Step1b_After_Petrogen
                            
                            RunConditions(runnumber,:)=StandardInfo;
                            
                         
                            
                            
                            
                            %  ForRevPet_Save(runnumber,:)=ForRevPet;
                            
                            runnumber=runnumber+1;
                        end
                        
                    end
                end
            end
        end
    end
end

runnumber = runnumber-1;







%%
% COMBINE the POOLS and the INSTANTS into 1 large matrix

% For each INFO matrix, make new matricies without the shared fun info (The
% first 30 columns), and make a new column that designates if the melt is a
% pooled melt or an incremental melt
STANDARD_INFO = MASTER_POOLS_INFO(:,1:size(StandardInfo,2));

MASTER_CAT_POOLS_INFO = MASTER_POOLS_INFO;
MASTER_CAT_POOLS_STDINFO=MASTER_CAT_POOLS_INFO(:,1:size(StandardInfo,2));
MASTER_CAT_POOLS_INFO(:,1:size(StandardInfo,2))=[];
MASTER_CAT_POOLS_type = [50015.*ones(size(MASTER_CAT_POOLS_INFO,1),1)]; %designates pooled melts


MASTER_CAT_column_POOLS_STDINFO=MASTER_column_POOLS_INFO(:,1:size(StandardInfo,2));
MASTER_CAT_column_POOLS_INFO=MASTER_column_POOLS_INFO(:,size(StandardInfo,2)+1:end);
MASTER_CAT_column_POOLS_type=[60150015.*ones(size(MASTER_CAT_column_POOLS_INFO,1),1)]; %designates pooled melts

MASTER_CAT_INSTANTS_INFO = MASTER_INSTANTS_INFO;
MASTER_CAT_INSTANTS_STDINFO(:,1:size(StandardInfo,2))=MASTER_CAT_INSTANTS_INFO(:,1:size(StandardInfo,2));
MASTER_CAT_INSTANTS_INFO(:,1:size(StandardInfo,2))=[];
MASTER_CAT_INSTANTS_type = [ 1115741175.*ones(size(MASTER_CAT_INSTANTS_INFO,1),1)];%designates instanteous melts

MASTER_CAT_SM_INFO = MASTER_SM_INFO;
MASTER_CAT_SM_STDINFO(:,1:size(StandardInfo,2))=MASTER_SM_INFO(:,1:size(StandardInfo,2));
MASTER_CAT_SM_INFO(:,1:size(StandardInfo,2))=[];
MASTER_CAT_SM_type = [size(PetrogenTraceElement_Strings,2).*ones(size(MASTER_CAT_SM_INFO,1),1)];%designates SM melts

% MASTER_SM_MAJORS_original=[MASTER_SM_MAJORS];
% MASTER_SM_TRACE_original=[MASTER_SM_TRACE];
% %MASTER_SM_INFO_original=[MASTER_SM_INFO];
% PadPoolsDs_SM_MAJORS =NaN.*zeros(size(MASTER_POOLS_MAJORS,1),size(MASTER_SM_MAJORS,2));
% PadPoolsDs_SM_TRACE =NaN.*zeros(size(MASTER_POOLS_MAJORS,1),size(MASTER_SM_TRACE,2));
% MASTER_SM_MAJORS =  [PadPoolsDs_SM_MAJORS; MASTER_SM_MAJORS];
% MASTER_SM_TRACE =[PadPoolsDs_SM_TRACE; MASTER_SM_TRACE];
% MASTER_SM_MAJORS
% MASTER_SM_TRACE
% MASTER_SM_INFO
% MASTER_column_POOLS_MAJORS
% MASTER_column_POOLS_TRACE


% Now pad a new matrix for each so that it has 17 columns, with pooled
% melts first,
%so the first column is the _type,
%then 10 entires of pooled info, then
%the 6 entires of instant info,
%and 4 entry of SM,
%and 4 entires of column pooled info
MASTER_CAT_POOLS_INFO = [MASTER_CAT_POOLS_STDINFO MASTER_CAT_POOLS_type...
    MASTER_CAT_POOLS_INFO ...
    NaN.*ones(size(MASTER_CAT_POOLS_INFO,1),6) ...
    NaN.*ones(size(MASTER_CAT_POOLS_INFO,1),4) ...
    NaN.*ones(size(MASTER_CAT_POOLS_INFO,1),4)];

MASTER_CAT_INSTANTS_INFO = [MASTER_CAT_INSTANTS_STDINFO MASTER_CAT_INSTANTS_type ...
    NaN.*ones(size(MASTER_CAT_INSTANTS_INFO,1),10) ...
    MASTER_CAT_INSTANTS_INFO ...
    NaN.*ones(size(MASTER_CAT_INSTANTS_INFO,1),4) ...
    NaN.*ones(size(MASTER_CAT_INSTANTS_INFO,1),4)];

MASTER_CAT_SM_INFO = [MASTER_CAT_SM_STDINFO MASTER_CAT_SM_type...
    NaN.*ones(size(MASTER_CAT_SM_INFO,1),10) ...
    NaN.*ones(size(MASTER_CAT_SM_INFO,1),6) ...
    MASTER_CAT_SM_INFO...
    NaN.*ones(size(MASTER_CAT_SM_INFO,1),4) ...
    ];

MASTER_CAT_column_POOLS_INFO=[MASTER_CAT_column_POOLS_STDINFO MASTER_CAT_column_POOLS_type...
    NaN.*ones(size(MASTER_CAT_column_POOLS_INFO,1),10) ...
    NaN.*ones(size(MASTER_CAT_column_POOLS_INFO,1),6) ...
    NaN.*ones(size(MASTER_CAT_column_POOLS_INFO,1),4) ...
    MASTER_CAT_column_POOLS_INFO];

MASTER_INFO = [MASTER_CAT_POOLS_INFO; MASTER_CAT_INSTANTS_INFO; MASTER_CAT_SM_INFO; MASTER_CAT_column_POOLS_INFO];
MASTER_MAJORS = [MASTER_POOLS_MAJORS; MASTER_INSTANTS_MAJORS; MASTER_SM_MAJORS; MASTER_column_POOLS_MAJORS];
MASTER_TRACE = [MASTER_POOLS_TRACE ; MASTER_INSTANTS_TRACE; MASTER_SM_TRACE; MASTER_column_POOLS_TRACE];

%MASTER_DISEQ = [MASTER_POOLS_DISEQ; MASTER_INSTANTS_DISEQ; NaN.*zeros(size(MASTER_SM_MAJORS,1),size(MASTER_POOLS_DISEQ,2))];
MASTER_columnPool_DISEQ = nan.*zeros(size(MASTER_column_POOLS_MAJORS,1),4);
MASTER_DISEQ = [MASTER_POOLS_DISEQ; MASTER_INSTANTS_DISEQ; MASTER_SM_DISEQ; MASTER_columnPool_DISEQ];


% keep original
MASTER_Residue_TRACE_original =[MASTER_Residue_TRACE];
MASTER_SolidResidue_TRACE_original =[MASTER_SolidResidue_TRACE];
MASTER_CPX_original =  [MASTER_CPX];
MASTER_INSTANTS_Ds_original =[MASTER_INSTANTS_Ds];

% make dummy matricies for pooled residues compositions (for now, could
%think about this more)
PadPoolsResidueTrace =NaN.*zeros(size(MASTER_POOLS_MAJORS,1)+size(MASTER_SM_MAJORS,1)+size(MASTER_column_POOLS_MAJORS,1),size(MASTER_Residue_TRACE,2));
PadPoolsSolidResidueTrace =NaN.*zeros(size(MASTER_POOLS_MAJORS,1)+size(MASTER_SM_MAJORS,1)+size(MASTER_column_POOLS_MAJORS,1),size(MASTER_SolidResidue_TRACE,2));
PadPoolsCPXTrace =NaN.*zeros(size(MASTER_POOLS_MAJORS,1)+size(MASTER_SM_MAJORS,1)+size(MASTER_column_POOLS_MAJORS,1),size(MASTER_CPX,2));
PadPoolsDs =NaN.*zeros(size(MASTER_POOLS_MAJORS,1)+size(MASTER_SM_MAJORS,1)+size(MASTER_column_POOLS_MAJORS,1),size(MASTER_INSTANTS_Ds,2));

%resave residue compositions
MASTER_Residue_TRACE =[PadPoolsResidueTrace; MASTER_Residue_TRACE];
MASTER_SolidResidue_TRACE =[PadPoolsSolidResidueTrace; MASTER_SolidResidue_TRACE];
MASTER_CPX =  [PadPoolsCPXTrace; MASTER_CPX];
MASTER_Ds =[PadPoolsDs; MASTER_INSTANTS_Ds];



%%
eval(sprintf('save(''%s'',''PetrogenBulkCompNames'',''PetrogenMajorElement_Strings'',''PetrogenTraceElement_Strings'',''Petrogen_Co'', ''Petrogen_Co_Name'', ''Petrogen_BulkCompositionsNormalized'', ''Petrogen_Dmatrix'',''MASTER_SM_DISEQ'',''MASTER_ALL_RESULTS'',''MASTER_SM_MAJORS'',''MASTER_SM_TRACE'',''MASTER_DISEQ'',''Petrogen_Dmatrix'',''outputMeltMatrix'',''outputGarnetMeltMatrix'',''outputMeltMatrix_FullPool'',''outputFirstMelt'',''MASTER_Ds'',''MASTER_INFO'',''MASTER_MAJORS'',''MASTER_TRACE'',''MASTER_SM_MAJORS'',''MASTER_SM_TRACE'',''MASTER_SM_INFO'',''MASTER_POOLS_INFO'', ''MASTER_POOLS_DISEQ'', ''MASTER_INSTANTS_DISEQ'',''MASTER_POOLS_MAJORS'', ''MASTER_POOLS_TRACE'', ''MASTER_INSTANTS_MAJORS'', ''MASTER_INSTANTS_TRACE'', ''MASTER_INSTANTS_INFO'',''MASTER_Residue_TRACE'',''MASTER_SolidResidue_TRACE'',''MASTER_CPX'',''MASTER_INSTANTS_Ds'')',[pwd '/' subfolder '/' worksheetName2Save '.mat']))



'Name of Folder Files are Saved in:'
subfolder
'Name of Saved Matlab Workspace and EXCEL spreadsheet:'
worksheetName2Save

%%
HeadingsPools = {'dFdPwaterchange','porosity','initialwatercontent','TP','U','BulkComp','FirstMeltPb','MeltExtractedPb','LastMeltPb','Gar2SpTransition','Sp2PlagTransition','aluminumOutPb','aluminumOutFrac','cpxOutPb','cpxOutFrac','Tdiff','crustalthickness','ct2_ininc(1)','percentGar','percentSpin','percentPlag','FP_Full','percent pool (1=full, 0 = on axis)',	'distance_ininc',	'Lnexp_cumulative_IT_ininc',	'%gar',	'%sp',	'%plag',	'gar km',	'sp km',	'plag km',	'totalkm'};


layers2plot=[1:31];

% high230Th = find(MASTER_POOLS_INFO(:,19)==0 & MASTER_POOLS_DISEQ(:,3)>1.2);
% mat2clip(MASTER_POOLS_INFO(high230Th,:))
%

toc


%%
if strcmp(IsobaricOrPolybaric,'Polybaric')
    disp('Making some useful P-T diagrams of melting...')
    if runnumber <= 4
        runnumberPTplot=runnumber;
    else
        runnumberPTplot=4;
    end
    
    axisdepth = 60;
    maxT = 1600;
    
    
    ConstraintOptions= {'dFdPwaterchange','rFrac_crit','Water','TP','U','Bulk Comp','FirstMeltPb','MeltExtractedPb','LastMeltPb','Gar2SpTransition','Sp2PlagTransition','aluminumOutPb','aluminumOutFrac','cpxOutPb','cpxOutFrac','Tdiff','quick crustal thickness','LastAdiabaticPressure','LastAdiabaticPressure2','LastAdiabaticPressure3','BrittleDuctile600Pres','ct2_ininc(1)','percentGar','percentSpin','percentPlag','percentGar_onaxis','percentSpin_onaxis','percentPlag_onaxis','FP_Full','FP_OnAxis','Fmax_OnAxis','Fmax_OnAxis_NotExp','230Th ZI FULL','226Ra ZI FULL','230Th CI FULL','226Ra CI FULL','maxNF 230Th ZI','maxNF 226Ra ZI','maxNF 230Th CI','maxNF 226Ra CI','TYPE','percentPool','distance_ininc','Variable Pool Fmean','Vpool %gar','Vpool  %sp','Vpool %plag','gar km','sp km','plag km','totalkm','Pb','T','Ts','Ln_cumulative','Lnexp_cumulative','StabilityField ','FN or FP','ZetaGar','ZetaSpinel','ZetaPlag','ColumnNum','variablypooledLnExp','variablypooledPb','variablypooledTs'};
    Variables4name = {'TP','Bulk Comp','Water','rFrac_crit','U'};
    [a,xxxvalue]=ismember(Variables4name,ConstraintOptions);
    
    for layerplotnum = 1:runnumberPTplot
        
        details = sprintf('%dC-#%d-%dppm-%s%d%%-%dcm/yr',...
            RunConditions(layerplotnum,xxxvalue(1)),RunConditions(layerplotnum,xxxvalue(2)),RunConditions(layerplotnum,xxxvalue(3)),'\phi',100*RunConditions(layerplotnum,xxxvalue(4)),RunConditions(layerplotnum,xxxvalue(5)));
        
        detailsv2 = sprintf('%dC-Bulk%d-%dppm-%dpercent-%dcm/yr',...
            RunConditions(layerplotnum,xxxvalue(1)),RunConditions(layerplotnum,xxxvalue(2)),RunConditions(layerplotnum,xxxvalue(3)),100*RunConditions(layerplotnum,xxxvalue(4)),RunConditions(layerplotnum,xxxvalue(5)));
        
        
        figure('units','normalized','outerposition',[0.0320 0.1777 0.2766 0.7285])
        FigureTitle = sprintf('Fig  P-T %s', detailsv2);
        set(gcf,'Name',FigureTitle)
        
        hold on
        
        layerplot_color_here = 'k';
        layerplot_marker = '*';
        plot(eval(sprintf('Ts_%g', layerplotnum)),Pbeek,'-','Color','k','linewidth',3)
        plot(eval(sprintf('T_%g', layerplotnum)),Pbeek,'-','Color','r','linewidth',3)
        %hline(eval(sprintf('FirstMeltPb_%g', layerplotnum)),'r')
        
        plot(eval(sprintf('Tb_%g', layerplotnum)),Pbeek,'-','Color','b','linewidth',3)
        plot(eval(sprintf('PooledT_%g', layerplotnum)),eval(sprintf('PooledPb_%g', layerplotnum)), 'ko-','MarkerSize',12);
        
        
        if any(eval(sprintf('FirstAt2percent_%g', layerplotnum)))>0
            yline(eval(sprintf('FirstAt2percent_%g', layerplotnum)),'r:','linewidth',3);
        end
        
        
        yline(eval(sprintf('LastMeltPb_%g', layerplotnum)),'r','Last Melt','LabelHorizontalAlignment','left','FontSize',12,'FontName','Times New Roman');
        
        Gar2SpTransition=Pb(eval(sprintf('Gar2SpTransition_index_%g', layerplotnum)));
        Sp2PlagTransition=Pb(eval(sprintf('Sp2PlagTransition_index_%g', layerplotnum)));
        
        
        yline(Gar2SpTransition,'k','Spinel','LabelHorizontalAlignment','left','FontSize',15,'FontName','Times New Roman');
        yline(Sp2PlagTransition,'k','Plagioclase','LabelHorizontalAlignment','left','FontSize',15,'FontName','Times New Roman');
        yline(Gar2SpTransition,'k','Garnet','LabelHorizontalAlignment','left','FontSize',15,'LabelVerticalAlignment','bottom','FontName','Times New Roman');
        
        
        
        
        yline(eval(sprintf('AveragePb_%g', layerplotnum)),'b','Average Pb on Axis','LabelHorizontalAlignment','left','FontSize',8,'LabelVerticalAlignment','bottom');
        
        yline(eval(sprintf('PooledPb_%g(1)', layerplotnum)),'b-','Full Pool','LabelHorizontalAlignment','left','FontSize',8,'LabelVerticalAlignment','bottom');
        
        
        %     yline(eval(sprintf('AllAveragePb_%g(1)', layerplotnum)),'c','Left Col','LabelHorizontalAlignment','left','FontSize',8,'LabelVerticalAlignment','bottom');
        %     yline(eval(sprintf('AllAveragePb_%g(end)', layerplotnum)),'c','Pool Right Col','LabelHorizontalAlignment','left','FontSize',12,'LabelVerticalAlignment','bottom');
        
        hleg = legend('Melting Curve','Mantle Temperature','Initial Mantle Temp','First extracted melt','location','best');
        title(hleg, details);
        
        
        axis([1100 maxT 0 axisdepth])
        ylabel('Pressure (kbar)')
        
        xlabel(sprintf('Temperature (%cC)',char(176)));
        box on
        %grid on
        set(gca, 'Xtick',200:100:1600)
        ax = gca;
        ax.XAxis.MinorTick = 'on';
        ax.XAxis.MinorTickValues = ax.XLim(1):50:ax.XLim(2);
        ax = gca;
        ax.YAxis.MinorTick = 'on';
        ax.YAxis.MinorTickValues = ax.YLim(1):1:ax.YLim(2);
        
        set(gca,'ticklength',1.5*[0.0200    0.0500])
        
        set(gca,'fontsize', 15,'LineWidth',0.7,'FontName','Times New Roman')
        set(gca,'XColor', 'k')
        set(gca,'YColor', 'k')
        
        
        set(gca,'YDir','Reverse')
        
        yyaxis right
        axis([1100 maxT 0 axisdepth*3.3])
        set(gca, 'Ytick',[0:10:axisdepth*33])
        ax.YAxis(2).MinorTick = 'on';
        ax.YAxis(2).MinorTickValues = [0:5:axisdepth*33];
        ax.YAxis(2).Color='k';
        ylabelax.YAxis(2).Color='k';
        set(gca,'YDir','Reverse')
        ylabel('Depth(km)')
        
    end
    
    % if JM==3
    % legend('Ts garnet','Ts H2O','Ts H2O dFdP','adiabat garnet','adiabat H2O','adiabat dFdP')
    % end
    %
    % if any(AluminumPhaseOut) ==1
    % hline(Pb(AluminumPhaseOut),'g')
    % text(1200, Pb(AluminumPhaseOut), 'Al phase out, dry')
    % end
    %
    % if any(AluminumPhaseOut_H2O) ==1
    % hline(Pb(AluminumPhaseOut_H2O),'g')
    % text(1200, Pb(AluminumPhaseOut_H2O), 'Al phase out, H2O')
    % end
    %
    % if any(AluminumPhaseOut_H2O_dFdP) ==1
    % hline(Pb(AluminumPhaseOut_H2O_dFdP),'g')
    % text(1200, Pb(AluminumPhaseOut_H2O_dFdP), 'Al phase out, H2O dFdP')
    % end
    %
    % if any(MeltExtractedPb_gar) ==1
    % hline(MeltExtractedPb_gar,'c-')
    % text(1200,MeltExtractedPb_gar,'critical melt fraction reached, dry')
    % hline(MeltExtractedPb_H2O,'c-')
    % text(1200,MeltExtractedPb_H2O,'critical melt fraction reached, H2O')
    % hline(MeltExtractedPb_H2O_dFdP,'c-')
    % text(1200,MeltExtractedPb_H2O_dFdP,'critical melt fraction reached, H2O dFdP')
    % end
    %
    
    
end


 Step1_Figures
end

